Dynamical analysis of the nonlinear growth of the m — n — 1 resistive internal mode 



CO 

o 
o 

CN 



M.-C. FirpcE and B. Coppi 
Massachusetts Institute of Technology, Cambridge, MA 02139-4307 
(Dated: February 2, 2008) 

A dynamical analysis is presented that self-consistently takes into account the motion of the 
critical layer, in which the magnetic field reconnects, to describe how the m = n = 1 resistive 
internal kink mode develops in the nonlinear regime. The amplitude threshold marking the onset of 
strong nonlinearities due to a balance between convective and mode coupling terms is identified. We 
predict quantitatively the early nonlinear growth rate of the m — n — 1 mode below this threshold. 
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The large scale dynamics and confinement properties 
of tokamak plasmas depend intimately on the behavior 
of m = n = 1 magnetohydrodynamic (MHD) internal 
kink modes. This has motivated an intense, long-lasting, 
experimental and theoretical research, notably devoted 
to study their implication in magnetic reconnection or as 
triggers of the sawtooth oscillations and crashes. These 
phenomena typically proceed beyond the linear regime, 
that is now rather well understood but assumes very 
small amplitudes of the modes. To offer a quantitative, 
predictive description of their nonlinear manifestations 
remains a difficult objective of both academic interest 
and very practical importance. This is especially relevant 
for the design of fusion burn experiments in which the ful- 
filment of linear stability constraints is challenged by the 
search for ignition. Such devices are thus expected to op- 
erate at best close to marginal stability for the m = n = 1 
ideal mode so that nonlinear effects come into play for 
fairly small values of the mode amplitude 0, 0] • 

In this Letter, we focus on the m = n = 1 resis- 
tive mode in which a finite resistivity ij destabilizes 
the otherwise marginally stable ideal MHD internal kink 
mode. Since Kadomtsev's scenario |^| predicting the 
complete reconnection of the helical flux within the q = 1 
surface on a timescale of order rf 1 !" 1 , that later appeared 
too large to account for observations, the nonlinear be- 
havior of the m = n = 1 mode has become a some- 
what controversial issue. Some numerical simulations 
suggested that the mode still grows exponentially into 
the nonlinear regime || which was supported by a the- 
oretical model p. Later some analytic studies 7\, sup- 
ported by numerical simulations rather predicted a 
transition to an algebraic growth early in the nonlinear 
stage. This result was challenged by Aydemir's recent 
simulations using a dynamical mesh [9|. These did show 
the linear exponential stage evolving towards an alge- 
braic stage, yet this was brutally interrupted by a second 
nonlinear exponential growth. A modified Sweet-Parker 
model was able to fit continuously both stages of evo- 
lution and the transition related to a change in the 
geometry of the current sheet [Hj- However, some funda- 
mental questions remain unanswered or unclear. Among 
them, how to relate the transition threshold with rj ? or 



what is the role of the (/-profile ? The aim of this Letter 
is to describe analytically how the m = n = 1 resistive 
mode develops in the nonlinear regime, by focusing on 
the equations controlling plasma dynamics. 

We consider the low-/? reduced MHD equations 



^ = [4>,U] + [JM 

^ = [<f>,1,]+T,(J-J ) 



(1) 

(2) 



assuming helical symmetry Only a single angu- 

lar variable is then involved in the problem, namely 
the helical angle a = ip — 0, with ip the toroidal 
and 9 the poloidal angles. U — Vj_<f> is the vortic- 
ity and J = Vj_^ the helical current density, with 
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d\. Time is normalized by the 



poloidal Alfven time (t — > t/rup), the radial variable 
by the minor radius (r — > r/a) and rj is the dimen- 
sionless resistivity, inverse of the magnetic Reynolds 
number S {rj = S^ 1 = thp/tr) with the poloidal 

1/2 

Alfven time th p — (moAj) R/Bq v and resistive time 
T R = MoQ 2 /vo- The Poisson brackets are defined by 
[4>, U] = -ip- (V_l</> x V ± U) = r- 1 (d r <j>d a U - d r Ud a cj>). 
4> and ip are the plasma velocity and helical magnetic 
field potentials expressed in cylindrical coordinates, so 
that the velocity is v = ip x V ±c/> and the magnetic field 
is B = B 0ip ip + 0xVx(tp~ r 2 /2). 

We consider MHD equilibria given by <po = and by 
an helical magnetic flux ipo (r) , related to the safety pro- 
file q(r) through d r ipo = r [1 — l/q(r)), such that q = 1 
for an internal radius r — r S Q. Thus d r ipo (fso) — 0. This 
means that the low-frequency ideal linear equations asso- 
ciated to are singular at r = r s o, with a formally 
diverging current density. This marks the presence of a 
critical layer in which the dynamics differs considerably 
from the outer one and where resistivity enters to cure 
the singularity. 

We wish to analyse perturbatively the time evolution 
of the m = 1 mode. For this, we assume that only the 
m = l mode is destabilized initially with an amplitude 
Aq, neglect all ideal MHD transients and restrict to the 
linear resistive timescale r = r\ x l z t. We do not consider 
the somehow ill-posed, singular limit r\ — > 0, but instead 
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realize that two small parameters are indeed competing 
in this problem, namely the small given resistivity 77 and 
the time-dependent amplitude A(t) of the linear m = 1 
mode. This introduces some subtleties in the amplitude 
expansion. The order one solution is given by linear the- 
ory using an asymptotic analysis |j| to match inner and 
outer solutions. Excitation of the m = 1 mode leads to 
a self-consistent correction to the location of the critical 
layer. One estimates the amplitude threshold, scaling 
with ?7, at which next order solution is required and the 
procedure iterated. Separability in time and space prop- 
agates at each order resulting in an amplitude expansion 
in A. As in any perturbative approach, the solution is 
formally known when the order one solution is. This is 
given by the linear theory reviewed now. 

Let /i, be the projection on exp(ima) of any function 
/at order n. In the inner resistive layer, Eqs. Q-0 read 



dr dx 2 1 
d 



■ d 2 /( D 



(3) 
(4) 



where we define k = ip' ' (r s o) /r s0 . In these equations, x 
is the stretched coordinate x = (r — r s0 ) /w and w = rj 1 / 3 
the magnitude of the width of the critical layer giving the 
maximal resistive ordering |3| in In the layer, 

radial derivatives are large, since d r — w~~ 1 d x and JSJ- 
Q are the dominant equations for w <C 1. There is one 
unstable solution, the m = 1 resistive mode, with growth 



rate 7l 



„ 2 /3 



l'(fso) 2 ^ 3 - Real space potentials read 



k 1/3 x 



ipi(x,a,T) = A exp(j L T) g L ^ j cos a (5) 

A Q ( k 1,3 x\ 

<pi(x,a,r) = --j=exp(j L T)g' L I -^=- j sina (6) 

where gi, is the function 

s 1 
9l (s) = - erfc(s)- ^=exp(-s 2 ). (7) 

This solution was chosen to satisfy the matching asymp- 
totic conditions g' L (— oo) = 1 and g' L (+oo) = 0. This 
analysis has to be complemented with the derivation 
of the outer solution. On the resistive timescale, this 
amounts to solve, at leading (zero) order in w, a lin- 
ear system of ideal MHD equilibria, singular at r — r s o 
[Tlj . This illustrates the passive character of the outer 
domain. We only retain here that, given the asymptotic 
and boundary conditions imposing ip± (rf ) = and 
ipi(l) = 0, the outer linear m = 1 solution ip^ '(r) is 
identically vanishing for r S Q < r < 1. 

Linear theory breaks down when, in the resistive crit- 
ical layer, nonlinear terms due to mode couplings, e.g. 



in Eq. Q [0i,i7i] ~ w 3 A 2 , balance linear ones, i.e. 
K xwd a Ji ~ A/w in Eq. ©. Thus A(t) = 0(r/ 2 / 3 ) 
marks the onset of second order terms. Before pursu- 
ing the analysis on the critical layer, we need to track 
it and self-consistently estimate its location. The total 
magnetic flux in the critical layer is now tp (x, a, t) = 
r] 2 / 3 ipQ (r s o) x 2 /2 + "01 { x i a i T )- To follow continuously 
the linear stage, we define the 'backbone' r s (a, r) of the 
critical layer as the 'neutral' field line with d r ip (r s ) = 0. 
Writing r s i(a, r) = r s (a,T) — r S Q = wxi{a,r) with 
d x tp(xi) = 0, this gives 



(a,r) 



Mr) 4 /3 gj(Q) 

V 1/3 yftVl (r s0 ) 



cos a 



(8) 



which relates to the shift of the core plasma inside the q = 
1 surface due to the kink instability. Then the x-point 
shift r s i (a = it, t) goes like Mt)/^ 1 / 3 , consistently with 
Aydemir's numerical results [9j. Thus the critical radius 
starts to leave the linear critical layer band, centered on 
r s o, when r s \ [a, r) becomes of the order r\ x l 3 for some 
a, that is when A (r) > rj 2 / 3 . This is again the threshold 
marking the end of the linear stage. We need now to 
define a generalized stretched coordinate in the critical 
layer as x = (r — r s (a, r)) jw. The replacements d T — > 



dr 



{dr s /di 



and d a — * d a — w 1 (dr s /da) d x 



are then required 

The second order critical layer equations involve an 
inhomogeneous part composed of quadratic terms in the 
order one solutions (J5J, © and iJSJ. This acts to force 
the growth of the m = and m = 2 perturbations but 
brings no contribution to the m — 1 dynamics. Therefore 
the m = 1 equations ©-0 are unchanged, except that, 
due to the motion of the critical layer JSJ, one needs to 
replace n n in ©-© by the time-dependent average 
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,(0) 



d 2 jj [r s (a, t)} 
(a,r) 



da. 



(9) 



This introduces a generalized linear system of equations. 
Neglecting the initially zero amplitudes of the m = and 
to = 2 perturbations in front of A(t), the second order 
correction to the location of the critical layer is given by 
r s2 (a, t) ~ - (2< (rso))- 1 C ( r *o) r.i (a, r) 2 . The va- 
lidity threshold of the second order solution is reached 
when the instantaneous critical line moves out of the 
critical layer of width w centered on r s0 + r sl (a, r) for 
some a. This corresponds to r s i (a,r) ~ w, that is to 
r sl (a,T) 2 ~ 77 1 / 3 , which gives A(t) = 0(r]^ 2 ). This 
threshold in the amplitude of the linear m = 1 mode 
marks the onset of third order terms, that will contribute 
again to the m = 1 dynamics. Its brutal manifestation 
is visible on Aydemir's plots 0- They clearly report a 
transition in the m = 1 kinetic energy when this becomes 
of order 77/2 [l^, namely around 5 x 1CP 8 for 7/ = 10~ 7 
and around 5 x 10~ 6 for 77 = 10~ 5 . 
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FIG. 1: Picture in the (x, a) space of the initial linear critical 
layer and of a nonlinear one centered on the instantaneous 
transverse neutral field line (in bold). The grey region repre- 
sents their overlapping domain within which the gradients of 
linear potentials are O(to- 1 )-large. 



At third order, cubic terms in the order one solutions 
or quadratic terms coupling the m — and m = 2 sec- 
ond order terms to the m = 1 first order ones appear 
in the inhomogeneous part of the critical layer equations 
and modify the m — 1 dynamics. These terms involve 
some radial derivatives, e.g. d r (j){ \ that are 0(w^ 1 )- 
large only within the linear layer. Locality enters here the 
analysis since the dominant contribution of these mode 
coupling terms comes from the localized zone in (r, a) 
where the instantaneous and linear critical layers over- 
lap. This is depicted by the grey shaded region in Fig. 
^ The novelty is that, in this region, mode couplings are 
now able to balance convective derivatives, both being 
dominant with respect to linear terms. More explicitly, 
while, e.g. in the Eq. written in the region where 
the instantaneous and linear critical layers overlap, the 
magnitude of linear terms is d T d 2 (f>^ ~ w~ 2 A(t), con- 
vective terms are of the order of d T r^ d^^-p ~ w~ 5 A 3 . 
Thus linear terms become negligible for A(t) ^> r? 1 ^ 2 , 
which marks the onset of the fully nonlinear regime for 
the m = 1 mode. Moreover, convective terms, e.g. 
^ • 1 - "~ 5 A 2 d T A, equilibrate mode coupling 
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terms, such as — d r <j)\~ 1S> d a r^ d r U^ 
ing from [0, U] in the shear- Alfven law (I}. The nonlinear 
growth rate derives from this balance. As k 1 -- ^ (t) is no 
longer involved in those convective and mode coupling 
terms, there is no extra time-dependence in the domi- 
nant equations, so that the nonlinear growth rate is just 
equal, by continuity, to the growth rate of the m = 1 
mode when A{r) becomes of order 77 1 / 2 . Its value de- 
pends notably on the equilibrium g-profile as we shall 
see below. After some spatial averaging, a rough sum- 
mary of the time evolution of the m = 1 mode amplitude 
may be then finally written as 



dA , . J c . 9 / dA . „ , 



= 0, (10) 



where the initial value of the growth rate 7 (0) is 7^ and 
where the early time dependence of 7 comes from the 
motion of the critical layer and is computed quantita- 
tively below. In Eq. (|10|) . c is a constant of order one 
and tffL denotes the (magnitude of the) time at which A 
becomes of order rj 1 / 2 . Eq. (fTU|) describes effectively the 
transition between two (almost) exponential stages. Be- 
cause 1^3 and ip^ are zero at the onset of the third order 
regime, Eq. (|1C)|> remains valid during some stage even if 
the structure and scaling of the critical layer should sub- 
stantially change as the generalized linear stage is left. 

For the convective exponential stage to be fully valid, 
the overlap between the linear and instantaneous critical 
layers should be large enough. One expects then a qual- 
itatively different late behavior of the m — 1 dynamics if 
the x-point region is far away from the linear layer when 
A(t) = O (V /2 ), that is, due to ©, if 77- 1 / 6 ^> 1. This 
regime is extremely challenging to reach numerically but 
may be satisfied in tokamak plasmas. 

We finally examine the early nonlinear effects on the 
growth rate of the m = 1 mode due to the motion of the 
critical layer. This amounts to solve the system of differ- 
ential equations for kq replaced with (r), de- 
fined in It can be checked that, as long as the order of 
magnitude of A (t) is lower than 77 1 / 2 , (r) may be ap- 
proximated by (2tt) 1 r s (a, r) 1 -0o [ r s (a, t)] da at 
leading order. This expression will be retained in the nu- 
merical computations. The time-dependent growth rate 
is defined as 7 (r) = d T A/A. In this generalized linear sys- 
tem, there is one condition shared with the linear deriva- 
tion: for a solution in separate variables t and x, it is 
that 7 (t) /k (t) be constant. This constant is then fixed 
by continuity with the linear solution at time zero giving 



7(t) 



7l 



-1/3 



(11) 



(r) 

Here one implicitly assumes that the spatial part of the 
linear eigenfunctions remains valid |14|. The instanta- 



neous critical radius is r s (a, r) 



f'sO 



V 



1/3, 



(a,- 



where x s (a,r) is given by the approximate expression 



x s (a,r) = H~ 



A(t)kI /3 



cos a 



772/372^(^0) 



(12) 



H 1 denotes the inverse of the monotonously growing 



k,q^ 3 x/V2 



Due to 



function defined by H(x) = x/g' L 
the asymmetric nature of the m = 1 resistive eigenfunc- 
tions 0, H~ 1 (x) is very asymmetric, grossly equal to 
x below x = and exponentially small above. This 
confers a much more important weight on negative 
arguments of H~ l than on positive ones in the averaging 
©. The magnetic island has thus a higher effective 
contribution to the early nonlinear correction of the 
growth rate than the region of x-point. A rough esti- 
mate of the angular average of x s is given by x^ 1 (r) ~ 
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FIG. 2: Analytic nonlinear growth rate corresponding to the 
initial conditions used in Ref. and resistivity r\ — 10 -7 , ne- 
glecting third order convective effects coming into play when 
A(t) becomes of order r/ 1 ^ 2 . This occurs for t ~ 1000. 
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FIG. 3: Analytic nonlinear growth rate for the same initial 
values as in Fig. |5] but with a modified equilibrium safety 
profile q(r). Its behavior around r s o is plotted in the insert. 



- (27T)- 1 A (r) kI' 3 / ( V ^V2^(r s0 )) dacosa. 
Eq. (|llfl defines a first order differential equation 
in A(t) that admits then the approximate form 

7 (r) ~ % + n^dr [r-^o W] foo) (r). Going back 
to time t and to = rf^lL, this gives 



C A(tf 



(13) 



where C = (q' + r sQ q'^ - 2r s0 qtf) / ^irV2r 2 s0 % /3 ^ and 

the index denotes an evaluation at r S Q. Eq. 1(13|) shows 
the first nonlinear contribution to the m = 1 evolution. 
The early behavior of the m = 1 growth rate is thus 
7(i) ~ 7l — CqAq exp ('jLt). In order to check numerically 
these analytic predictions for the generalized linear stage, 
that brings the first nonlinear contributions to the growth 
rate, we used Aydemir's initial conditions 0. The safety 

r a r 2 1 1 1/2 

profile is q(r) = q m jl + r 4 (q a /q m ) _1 J with 
g m = 0.9, g a = 3, giving Co > 0. The differential equa- 
tion l |l 1(1 was integrated numerically for Aq = ^/2x 10~ 55 
corresponding to an initial kinetic energy in the m = 1 
mode of the order 10 -11 . The nonlinear growth rate 
7(i) = Ty 1 / 3 ^) is plotted on Fig. Hfor S = 10 7 . This 
curve appears to be in fine agreement with the Figure 1 
of Ref. 9] for times t roughly below 1000 Alfven times. 

Fig. |31 illustrates the influence of the g-profile around 
r s0 on the time evolution of 7 due to © . A sudden bump 
in the nonlinear growth could thus even be observed, be- 
fore the onset of convective effects, for the special shape 
of q chosen in Fig. |3| Moreover, some g-profile may in- 
duce a saturation of A below the convective threshold and 
lead to partial reconnection. Most importantly, the ap- 
proach described here may be transposed to model the 
early nonlinear behavior of a variety of internal kinks 
such as two- fluid 0, 0, 0] and/or collisionless 
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